﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace RaytracerLib.MathLib.Polynomials
{
    [DebuggerStepThrough]
    public struct SturmSequence
    {
        private readonly Polynomial[] m_seq;
        public readonly int Length;

        public SturmSequence(Polynomial a_poly)
        {
            a_poly = a_poly.RemoveZeroes();

            m_seq = new Polynomial[a_poly.Order + 1];

            Length = 1;
            m_seq[0] = a_poly;

            if (a_poly.Order == 0)
                return;

            Length++;
            m_seq[1] = a_poly.Differentiate();

            while (m_seq[Length - 1].Order != 0)
            {
                Length++;
                m_seq[Length - 1] = -(m_seq[Length - 3] % m_seq[Length - 2]);
                if (m_seq[Length - 1].Order == 0)
                    return;
            }
        }

        public static Polynomial SingleRoots(Polynomial a_poly, double a_precision = Constants.DOUBLE_PRECISION)
        {
            return a_poly / Polynomial.GCD(a_poly, a_poly.Differentiate(), a_precision);
        }

        public Polynomial this[int a_index]
        {
            get
            {
                if (a_index >= Length)
                    throw new IndexOutOfRangeException();

                return m_seq[a_index];
            }
        }

        public int NumberOfSignChangesAt(double a_x)
        {
            int sc = 0;

            double c1 = m_seq[0].Evaluate(a_x);

            for (int i = 1; i < Length; i++)
            {
                double c2 = m_seq[i].Evaluate(a_x);

                if ((c1 * c2 < 0) || (c1 == 0.0))
                    sc++;

                c1 = c2;
            }

            return sc;
        }

        public int NumberOfSignChangesAtPositiveInfinity()
        {
            int sc = 0;

            double c1 = m_seq[0].LeadCoefficient;

            for (int i = 1; i < Length; i++)
            {
                double c2 = m_seq[i].LeadCoefficient;

                if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
                    sc++;

                c1 = c2;
            }

            return sc;
        }

        public int NumberOfSignChangesAtNegativeInfinity()
        {
            int sc = 0;

            double c1 = -m_seq[0].LeadCoefficient;
            if (m_seq[0].Order % 2 == 0)
                c1 = -c1;

            for (int i = 1; i < Length; i++)
            {
                double c2 = -m_seq[i].LeadCoefficient;
                if (m_seq[i].Order % 2 == 0)
                    c2 = -c2;

                if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
                    sc++;

                c1 = c2;
            }

            return sc;
        }
    }
}
